%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nround = 1e4;

z1_coint1pos = zeros(lent - 1, nround);
z1_coint2pos = zeros(lent - 1, nround);
z1_spos = zeros(lent - 1, nround);

parfor iround = 1:nround
    X1 = zeros(size(X2));
    shock_rng = randi([1, T - 1], 1, T);

    for t = 2:T
        X1(t, :) = X1(t - 1, :) * Psi' + eps(shock_rng(t), :);
    end

    [Psi1, Sig1] = run_VAR(X1, N, T, inflpos, y1pos, ysprpos, gdppos, divgrpos, pdpos, dtpos, dgpos, ...
        coint_tax, coint_spending);

    z1 = zeros(lent, N);
    z1(2, :) = Sig(:, shocks(1))' / Sig(shocks(1), shocks(1)) * shock_size(1);

    for t = 3:lent
        z1(t, :) = Psi1 * z1(t - 1, :)';
    end

    spos = N + 1;
    z1(:, spos) = mean(taxrevgdp) * exp(z1(:, coint_tax)) - mean(spendgdp) * exp(z1(:, coint_spending));

    z1_coint1pos(:, iround) = (mean(taxrevgdp) * (exp(z1(2:lent, coint_tax)) - 1)) * 100;
    z1_coint2pos(:, iround) = (mean(spendgdp) * (exp(z1(2:lent, coint_spending)) - 1)) * 100;
    z1_spos(:, iround) = (z1(2:lent, spos) - z1(1, spos)) * 100;

    z2 = zeros(lent, N);
    z2(2, :) = Sig(:, shocks(2))' / Sig(shocks(2), shocks(2)) * shock_size(2);

    for t = 3:lent
        z2(t, :) = Psi1 * z2(t - 1, :)';
    end

    spos = N + 1;
    z2(:, spos) = mean(taxrevgdp) * exp(z2(:, coint_tax)) - mean(spendgdp) * exp(z2(:, coint_spending));

    z2_coint1pos(:, iround) = (mean(taxrevgdp) * (exp(z2(2:lent, coint_tax)) - 1)) * 100;
    z2_coint2pos(:, iround) = (mean(spendgdp) * (exp(z2(2:lent, coint_spending)) - 1)) * 100;
    z2_spos(:, iround) = (z2(2:lent, spos) - z2(1, spos)) * 100;

    z3 = zeros(lent, N);
    z3(2, :) = Sig(:, shocks(3))' / Sig(shocks(3), shocks(3)) * shock_size(3);

    for t = 3:lent
        z3(t, :) = Psi1 * z3(t - 1, :)';
    end

    spos = N + 1;
    z3(:, spos) = mean(taxrevgdp) * exp(z3(:, coint_tax)) - mean(spendgdp) * exp(z3(:, coint_spending));

    z3_coint1pos(:, iround) = (mean(taxrevgdp) * (exp(z3(2:lent, coint_tax)) - 1)) * 100;
    z3_coint2pos(:, iround) = (mean(spendgdp) * (exp(z3(2:lent, coint_spending)) - 1)) * 100;
    z3_spos(:, iround) = (z3(2:lent, spos) - z3(1, spos)) * 100;
end
